Bouncing ball orbits and symmetry breaking effects in a three-dimensional chaotic 

billiard 
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We study the classical and quantum mechanics of a three-dimensional stadium billiard. It consists 
of two quarter cylinders that are rotated with respect to each other by 90 degrees, and it is classically 
chaotic. The billiard exhibits only a few families of nongeneric periodic orbits. We introduce an 
analytic method for their treatment. The length spectrum can be understood in terms of the 
nongeneric and unstable periodic orbits. For unequal radii of the quarter cylinders the level statistics 
agree well with predictions from random matrix theory. For equal radii the billiard exhibits an 
additional symmetry. We investigated the effects of symmetry breaking on spectral properties. 
Moreover, for equal radii, we observe a small deviation of the level statistics from random matrix 
theory. This led to the discovery of stable and marginally stable orbits, which are absent for unequal 
radii. 

PACS numbers: 05.45.-a, 03.65.Sq, 41.20.-q, 41.20.Jb 



I. INTRODUCTION 

Wave chaotic phenomena are visible in a large variety 
of physical systems ranging from lattice QCD [l], [2] to 
nuclei Q, atoms [J], mesoscopic systemslaL optical mi- 
crocavities @ , microwave resonators 0, H, Q , and to vi- 
brations of macroscopic objects [Ty, [Tl|]. The correspon- 
dence between the classical dynamics and wave phenom- 
ena in the semiclassical regime is of particular interest in 
such systems; for a comprehensive review, we refer the 
reader to Ref. [l2|, Q3 • It is best understood in Hamilto- 
nian systems with two degrees of freedom, whereas there 
is a lack of studies of chaotic three-dimensional systems. 
Experimental investigations of wave chaotic phenomena 
in three dimensions can be performed with microwave 
resonators (13, lla. Il6] and acoustic blocks (HALLO]- The 
resonance spectra investigated in such experiments are 
described by non-scalar wave equations and are consid- 
erably more complicated than the quantum mechanical 
wave equations of a non-integrable Hamiltonian system 
with three degrees of freedom. 

Let us briefly review the (relatively short) list of studies 
of quantum chaos in three-dimensional billiards. Prosen 
investigated quantum chaotic phenomena in a three- 
dimensional deformed sphere [17] . Primack and Smilan- 
sky [18] unraveled the classical and quantum mechanics 
of the three-dimensional Sinai billiard [19] . and verified 
the applicability of Gutzwiller's trace formula [20|, 121] : 
a corresponding experimental study of a microwave res- 
onator was presented in Ref. [16] , The experimental 
study J22j | of the three-dimensional Bunimovich stadium 
led to a verification of the trace formula proposed by 
Balian and Duplantier [23] for electromagnetic systems. 
A self-bound three-body system with high-dimensional 
scars was studied theoretically in Ref. [24] . Problems in- 



volving mode coupling in three-dimensional systems were 
investigated in vibrating crystals jll] and also in a mi- 
crowave resonator [22] . 

The theoretical description of quantum chaos in three- 
dimensional systems is in general very difficult. A full 
understanding of the classical dynamics is compounded 
by the difficulty to visualize motion in phase space. More- 
over, the analysis of level statistics of the corresponding 
quantum system requires the accurate computation of 
long level sequences and can be a very time-consuming 
task for eigenstates with short wave lengths. Further- 
more, there are a number of open questions. These con- 
cern, for instance, the applicability and accuracy of semi- 
classical periodic orbit sums for the quantization of such 
systems. This problem was addressed by Primack and 
Smilansky in their study of the three-dimensional Sinai 
billiard [l8|, which is completely chaotic. Yet, an infinite 
number of families of marginally stable (bouncing ball) 
orbits considerably complicates the semiclassical compu- 
tation of the level density in terms of Gutzwiller's trace 
formula. For the evaluation of the resonance density 
in three-dimensional microwave resonators, Gutzwiller's 
trace formula does not apply, and one has instead to use 
the trace formula by Balian and Duplantier. The ap- 
plicability of this periodic orbit sum was recently inves- 
tigated and confirmed for a three-dimensional stadium 
billiard with chaotic dynamics (22(]. It is the purpose of 
the present work to study the quantum mechanical as- 
pects of the three-dimensional stadium billiard further. 
We will focus in particular on quantum manifestations of 
classical chaos, the applicability of Gutzwiller's periodic 
orbit sum, and on effects related to symmetry breaking. 

The billiard depicted in Fig. [1] is a generalization of 
the two-dimensional Bunimovich stadium [25] to three 
dimensions [26| . It consists of two quarter cylinders that 



are rotated about 90 degrees with respect to each other. 
For the case that these quarter cylinders are separated 
by a finite distance a, Bunimovich and Del Magno |27j 
showed that this billiard is completely hyperbolic, i.e. 
there is no finite measure of trajectories that is not ex- 
ponentially sensitive to changes of their initial condi- 
tions. Earlier numerical studies suggest that the billiard 
of Fig. Q]is also classically chaotic [26j. In what follows, 
we restrict ourselves to this case. 




FIG. 1: (Color online) Three-dimensional stadium billiard. 
The plane z — is shaded. 

There are several reasons why the three-dimensional 
stadium billiard is a promising candidate for the study of 
quantum chaos. First, the classical dynamics is strongly 
chaotic and no stable islands are known for v\ ^ r 2 . In 
this work, we focus particularly on the case n = v2^2 
since this geometry was studied in microwave experi- 
ments |22| and in Ref. [2y|. The ratio was chosen irra- 
tional in order to avoid nongeneric quantum effects due 
to classical orbits of measure zero (see below). Second, 
nongeneric contributions to the level density arise due to 
two families of bouncing-ball orbits. However, unlike in 
the case of the three-dimensional Sinai billiard [18} , their 
contribution can be computed and subtracted. Third, 
for r\ — r 2 the billiard exhibits a high symmetry (under 
reflection at the z = plane with a subsequent rotation 
around the z-axis about 90 degrees), and this offers the 
opportunity to study symmetry-breaking effects. Below, 
we will see that the level statistics in this highly sym- 
metric case exhibit a deviation from the theoretical pre- 
diction for chaotic systems. This led to the discovery of 
a stable and a few marginally stable orbits which were 
previously unknown. 

This article is organized as follows. In Sect. HH wc 
focus on classical periodic orbits and nongeneric modes. 
In Sect. IIII1 we present our analysis of the quantum me- 
chanical spectra and discuss level statistics and symme- 
try breaking. We conclude with a summary of the main 
results. Some technical aspects concerning the finite ele- 



ment approximation with web-splines are deferred to the 
Appendix. 



II. NONGENERIC MODES AND PERIODIC 
ORBITS 

In this section we investigate nongeneric modes and 
classical periodic orbits. The former are an interesting, 
system-specific property and must be subtracted from 
the staircase function before generic properties can be 
analyzed; the latter are useful in a semiclassical interpre- 
tation of quantum spectra. 



A. Nongeneric modes 

We study the desymmetrized stadium billiard depicted 
in Fig. [1] The dynamics is limited to x > and y > 
with specular reflections at the planes x = and y = 0. 
We use ?"i > r 2 and will express lengths in units of r 2 and 
wave momenta in units of r^ . For r\ ^ r 2 , the billiard 
is fully desymmetrized; for r\ = r%, it is still symmet- 
ric under reflection at the plane z = and a subsequent 
rotation by ir/2 around the z-axis. We label eigenstates 
by their wave momentum k. Employing a recently de- 
veloped method which is based on finite elements and 
seeks solutions of the Schrodinger equation with Dirich- 
let boundary conditions, we computed the lowest 1200 
levels up to fcr 2 ~ 35. More details on this method are 
given in the appendix. 

The staircase function 



x 



N(k) = ^2@(k-k n ) 



(1) 



with Q(x) denoting the unit step function counts the 
levels below a given k. It is the sum of a smooth 
part N smoot h(k) and a fluctuating part iVa uc (fc), i.e., 
N(k) — N smoo th(k) + Nf{ uc (k). The smooth part is given 
by the Weyl formula which is a polynomial of degree three 
in k, and is subtracted from N(k). Figure [2] shows the 
remaining fluctuating part. One identifies large-scale os- 
cillations that grow in amplitude with increasing wave 
momentum. These fluctuations are due to the nongeneric 
modes of the billiard associated with the bouncing-ball 
orbits perpendicular to the flat boundaries of the billiard 
and the orbits within the z = plane. 

The fluctuating part of the level density is obtained by 
differentiation, i.e. ptiuc(k) = ^riVfl uc (A;). The absolute 
value squared of the Fourier transform of this quantity 
yields the length spectrum. It is shown in the upper 
part of Fig. [3] for n = v / 2r 2 . The peaks in the length 
spectrum appear at the lengths of classical periodic or- 
bits [2l[. The most prominent peaks are at lengths that 
are multiples of the bouncing-ball orbits with lengths 1t\ 
and 2r 2 , respectively. There are other peaks that can be 




FIG. 2: (Color online) Fluctuating part iVfl uc of the staircase 
function for the geometry with ri/ra = v2 (thin line). Con- 
tributions from bouncing-ball orbits (dashed line) and from 
all nongeneric modes (thick line). 



identified with the lengths of orbits inside the rectangu- 
lar z — plane. Both types of orbits are nongeneric 
due to their particular stability properties. As we are 
interested in generic properties, we have to extract these 
contributions. 







FIG. 3: (Color online) Upper part: Length spectrum of the 
three-dimensional stadium billiard for n/ra = v2 (full line). 
Contributions from unstable periodic orbits (dashed line). 
Lower part: Sum of contributions from nongeneric modes and 
unstable periodic orbits (full line). Peaks at multiples of 2n 
and 2ri correspond to the bouncing-ball orbits in the two 
quarter cylinders. Diamonds: Lengths of nongeneric orbits in 
the 2 = plane. 

First, let us consider the bouncing-ball orbits. There 
are two families of bouncing-ball orbits which bounce 
back and forth between the two parallel planes of dis- 
tance n and r2 in the two quarter cylinders, respec- 
tively. These orbits are marginally stable as they have 
zero Lyapunov exponents. According to semiclassical pe- 



riodic orbit theory, they are associated with long-range 
oscillations of the quantum mechanical density of states. 
We follow Ref . [15| for the computation of the contribu- 
tion of bouncing-ball orbits to the staircase function, and 
focus on the bouncing-ball modes in the upper quarter 
cylinder first. The number of bouncing-ball (bb) modes 
up to momentum k is given by 
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and l x = ri, l y = r%. This ansatz for the number 
of bouncing-ball modes is intuitively clear. The trace 
is performed in the semiclassical approximation as a 
phase space integral in the directions perpendicular to 
the bouncing-ball orbits, while these modes are quantized 
along the orbit. The integrations can be performed, and 
one obtains [l5[ 



iVbbW - ^^(fc 2 -fcye(fc 2 -fc : 
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^E( fc2 - fc U ( fc2 - fc U • ( 4 ) 



To obtain the fluctuating contribution iVbb,fluc, we sub- 
tract a polynomial of third degree in k from the expres- 
sion (0| . We confirmed that the leading coefficients of 
this polynomial agree with the coefficients of the Weyl 
formula. An analytical formula is presented in Ref. [1 51 ] - 
As evident from Fig. [51 the bouncing-ball orbits (dashed 
line) describe the gross oscillations of the fluctuating 
staircase function rather well. To account for finer de- 
tails, we also have to consider nongeneric modes inside 
the z — plane. 

The rectangular z = plane is special, since the classi- 
cal motion is regular inside this plane but unstable with 
respect to deviations out of this plane. Note that this 
does not spoil the hyperbolicity of the billiard since the 
set of orbits restricted to the z = plane are of mea- 
sure zero in classical phase space. However, remnants of 
families of classical orbits inside this plane can also be 
associated with long-range fluctuations in the quantum 
mechanical level density for the following two reasons. 
First, the periodic orbits inside this plane come in two- 
dimensional families opposed to the isolated periodic or- 
bits in hyperbolic systems. Second, these orbits are less 
unstable than most truly three-dimensional periodic or- 
bits and thus enter any periodic-orbit sum with a greater 



weight. It is thus necessary to include the nongeneric 
modes associated with the z = plane. 

Extending the results from Ref. [l5j, the number of 
modes associated with nongeneric orbits is given by the 
quantum mechanical trace over nongeneric (ng) modes 



N ng (k) = Tr ng 8 



dzdk z 
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This ansatz is again transparent: the trace is approxi- 
mated by a phase-space integral in the ^-direction, while 
the rectangular sections (with z-dependent area) perpen- 
dicular to the z-axis are explicitly quantized. We thus 
quantize the motion inside this plane adiabatically, quite 
similar to the approach that Bai et al. [28fl applied to the 
Bunimovich stadium. Note that Eq. (0 accounts both 
for the bouncing-ball orbits and for the modes associated 
with the z = plane. 

In Eq. (J5|), the wave momenta k x .^ and ky v are defined 



as in Eq. ([3]) but with the z-dependent radii l x 
We have 



and 



lx(z) = 



l y (z) = 



r\ for z < 0, 

yjr\ — z 2 for z > 0. 

r 2 for z > 0, 

\Jr\ — z 2 for z < 0. 



(6) 



The integrations in Eq. ([5]) can be performed analytically, 
and one obtains 

N ng (k) = -Y] {A(K^,r 1 ,r 2 ) + A(K IM/ ,r 2 ,r 1 )) . (7) 

7T ' ' 
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Here 



A(Kpv,ri,r 2 ) = r x ^j k^ - k\ v 

(E( K/1V ) - (1 - /&,)#(«„„)) (8) 



with 

K^v = K IMy (k,r 1 ,r 2 ) 



k 2 -k 2 J 



(9) 



contains the dependence on the wave momentum and the 
radii of the billiard. The wave momenta k\ t ^ and k 2%v 
are defined in Eq. ([3]). In Eq. ||5J), K (k) and E(k) are the 
complete elliptical integrals of the first and second kind, 
respectively. The prime in the sum of Eq. ([7]) indicates 
that the summation is limited to values of /z, v such that 
the square roots in Eqs. |(5J| and ([9]) are real. 

The smooth part of the staircase function is again a 
polynomial of degree three and is subtracted. The fluc- 
tuating part of N ng (k) is shown in Fig. [5] as a thick line. 
Clearly, the gross and also finer oscillations are accounted 



for. This enables us to analyze the length spectra de- 
picted in Fig. [3J The upper part shows the length spec- 
tra, i.e. the Fourier transform of the fluctuating part 
of the density of states, the lower part that of the non- 
generic modes and the unstable periodic orbits, which 
has been computed based on Eq. ((7j) and the Gutzwillcr 
trace formula (see next section). The peaks at multiples 
of the lengths 2r\ and 2r 2 are due to the bouncing-ball 
modes. The peaks marked by a diamond are due to the 
nongeneric orbits in the z = plane. The remaining 
peaks must thus be associated with lengths of isolated 
periodic orbits, and we turn to their analysis in the fol- 
lowing subsection. 



B. Unstable periodic orbits 

Following Gutzwiller's periodic orbit theory [2l| . the 
semiclassical approximation of the quantum mechanical 
density of states is given in terms of the periodic orbits of 
the underlying classical system. This sum is infinite such 
that approximations have to be invoked. Here, we are 
interested in the length spectrum, i.e. the power spec- 
trum of the Fourier transform of the fluctuating part of 
the level density. The length spectrum up to length I is 
given in terms of all periodic orbits up to length I, and 
this is a finite (but usually with I exponentially increas- 
ing) number of periodic orbits. 

The search for periodic orbits is a cumbersome task. 
Here, we use two different methods and focus on peri- 
odic orbits outside the z = plane. The first method 
considers the Poincare surface of section (PSOS) defined 
by z — 0,p z > 0, and constructs the PSOS map. Peri- 
odic orbits are fixed points of this map. They are found 
by starting a large number of trajectories in the PSOS 
and by using a Newton-Raphson algorithm to find fixed 
points in the vicinity of each such trajectory. 

The second method is based on a symbolic code and 
utilizes the fact that the length of a periodic orbit is a 
local extremum under variation of the points of reflection 
along the billiard's boundary. This procedure is similar 
to the one described in Ref. [18j. We consider an open 
billiard system which is the infinite periodic extension 
of our billiard, and assign the letters "+" and "-" to 
reflections on the curved parts of the upper and lower 
quarter cylinder, respectively. Periodic orbits outside the 
z = plane can certainly be described as words composed 
from these two letters. We have no proof that there is 
a one-to-one correspondence between this symbolic code 
and the periodic orbits of our system. Given a word 
from this alphabet, we construct a random closed orbit as 
follows. For each "+" ("-") letter of the word, we choose 
a random point on the curved surface of the upper (lower) 
quarter cylinder. We connect this sequence of points by 
straight lines and compute the length of this closed orbit. 
Then we vary the positions of the random points in order 
to find a local minimum of the length. 

Both methods yield a considerable number of periodic 



A+A_ = 
with A + A 



orbits, and we consider the union of the resulting sets 
as our set of periodic orbits. Once a periodic orbit is 
found, we compute its length, monodromy matrix and 
its Maslov index following Ref. [29|. Recall that the 
monodromy matrix is the tangent map and encodes the 
stability properties of the periodic orbit. For unstable 
orbits, its eigenvalues come in real pairs A + ,A_ with 
1 or in complex quadruples A + , A_, A5_, A* 
_ = 1. Stable orbits have two pairs of com- 
plex eigenvalues with modulus one. 

We performed our most extensive search for the bil- 
liard with the geometry n = V2r 2 and determined more 
than 2000 periodic orbits. We believe that our list is 
fairly complete for the shortest orbits. All periodic obits 
we found were unstable. This extends and confirms the 
numerical results of Ref. [26| and strongly suggests that 
the billiard is completely chaotic. The quantum mechan- 
ical results presented below support this picture. The 
contribution of the unstable periodic orbits to the length 
spectrum is shown in the upper part of Fig.[3]as a dashed 
line. Clearly, unstable periodic orbits contribute signif- 
icantly to the length spectrum, particularly for length 
Ijr-i > 6 where the contributions from nongeneric modes 
become less dominant. 

However, for the case v\ = r 2 we found a stable peri- 
odic orbit. This was unexpected, since earlier studies [26| 
yielded a positive Lyapunov exponent for a sample of 10 4 
randomly chosen trajectories, and not one of the trajecto- 
ries was found to be stable. The shortest stable periodic 
orbit we found has a length l/r 2 ~ 10.1706. Note that 
the eigenvalues of the monodromy matrix corresponding 
to this orbit are on the unit circle and very close to the 
real axis. Thus it is difficult to distinguish this orbit from 
a marginally stable orbit. We also obtained periodic or- 
bits (in addition to the bouncing-ball orbits) that are 
marginally stable, i.e. the eigenvalues of the monodromy 
matrix are real with modulus one. The shortest of these 
orbits has a length l/r 2 ~ 5.4981. Other marginally sta- 
ble orbits have periods l/r 2 ~ 5.5153,6.2701,6.4738. 

To study the relevance of these orbits, we show the 
length spectrum of the billiard for n = r 2 in Fig. 2] Di- 
amonds denote the lengths of bouncing-ball orbits and 
of nongeneric periodic orbits associated with the z — 
plane; they account for many peaks in the length spec- 
trum. We also computed the lengths, monodromy ma- 
trices and Maslov indices of the 400 shortest unstable 
and marginally stable orbits by proceeding as in the case 
n = V2r 2 . However, we have no semiclassical the- 
ory for the computation of the density of states for the 
marginally stable and the stable orbit, and we are not 
able to disentangle their contribution to the length spec- 
trum from that of the nongeneric modes. For this reason 
we show in Fig. 3] the length spectrum of the billiard 
and indicate the lengths of the nongeneric orbits by dia- 
monds, and those of the marginally stable orbits by the 
full arrows. The dashed arrow marks the length of the 
stable orbit. It is difficult to clearly identify its impact 
on the length spectrum, since the peak around length 



l/r 2 m 10.19 can be attributed to the stable orbit and/or 
to a nongeneric orbit. Recall that a stable orbit leaves 
a strong imprint in the length spectrum if the stable is- 
land around it is sufficiently large (in units of (27rfi,) 3 ) 
to accommodate eigenstates. To obtain an estimate for 
the phase-space volume of the elliptical island we started 
bundles of 21 4 trajectories in the PSOS close to the sta- 
ble periodic orbit. All trajectories departed from the sta- 
ble orbit after a few intersections with the PSOS. Within 
the achievable numerical accuracy we may thus conclude, 
that either the volume of the phase space associated with 
the stable orbit is very small or that contrary to our nu- 
merical results it is only marginally stable. Note that 
there are several marginally stable orbits associated with 
the third arrow from the left. These orbits have almost 
identical lengths, and the visual inspection indicates that 
they are whispering gallery orbits. In this case, interfer- 
ence effects might explain why there is no clear peak in 
the length spectrum associated with these orbits. A sim- 
ilar effect has been found in the two-dimensional stadium 
billiard [3(| ■ In the next section we present evidence that 
the stable and marginally stable orbits explain peculiar- 
ities in the level statistics of the billiard with n = r 2 . 




FIG. 4: (Color online) Full line: Length spectrum for n — t%. 
Diamonds: Length of bouncing-ball orbits and of nongeneric 
orbits in the z = plane. Full arrows: Length of marginally 
stable orbits. Dashed arrow: length of stable periodic orbit. 



III. SPECTRAL PROPERTIES AND 
SYMMETRY BREAKING 

In this section we investigate statistical properties of 
the eigenvalues of the three-dimensional stadium billiard 
for a varying ratio of the radii r\/r 2 of the two quar- 
ter cylinders, where the volume of the billiard is kept 
fixed. We computed the lowest 1200 levels with a re- 
cently developed finite element method [3l| described in 
the appendix. The lowest 200 levels are discarded since 



we are interested in the correspondence between classical 
and quantum mechanics. In order to get rid of system- 
dependent properties, we need to rescale the eigenvalues 
to unit mean spacing and to extract the contributions 
of the nongeneric modes to the fluctuating part of the 
resonance density. This is done by replacing the com- 
puted wave numbers k by k = N sraoo ^(k) + A r lis .n, lc (/r) 
(see [9, [l5J]). As the underlying classical dynamics is 
completely chaotic, we expect agreement of the statisti- 
cal properties of the unfolded eigenvalues with those of 
random matrices drawn from the Gaussian orthogonal 
ensemble (GOE) if r\ is chosen not equal to r 2 32]. In 
Fig. [5]we show the nearest- neighbor spacing distribution 
and the S 2 -statistics for a ratio of radii r\jr% = \/2, 
i.e. for the case considered in the experiments with the 
microwave cavity of the shape of a three-dimensional 
stadium billiard [22 |. Both curves agree well with the 
corresponding ones for random matrices from the GOE. 
This shows that first the three-dimensional quantum sta- 
dium billiard behaves like a generic quantum system with 
chaotic classical dynamics, and second that our proce- 
dure of extracting the contribution of nongeneric modes 
is applicable and complete. 

We recall that the experiment [221 ] with a microwave 
resonator of the shape of a three-dimensional stadium 
billiard with n /r% = \/2 revealed deviations of the spec- 
tral properties from GOE-behavior. This was attributed 
to a partial decoupling between electromagnetic TE and 
TM modes in the low-frequency regime. Note, however, 
that this system differs from the quantum billiard con- 
sidered in the present work due to the vectorial character 
of the underlying wave equations. As a consequence, the 
trace formula by Balian and Duplantier [23] has to be 
employed for the semiclassical computation of the length 
spectrum. It differs from Gutzwiller's trace formula in 
the occurrence of a factor associated with the polariza- 
tion, and as a consequence, only periodic orbits with an 
even number of reflections enter. The evaluation of this 
factor showed that the electric and magnetic modes are 
decoupled on a considerable fraction of the short peri- 
odic orbits. This finding supports the interpretation of 
the observed deviations from GOE-behavior in terms of 
a partial decoupling between the TE and the TM modes. 
In [221 ] the experimental length spectrum was well repro- 
duced by the theoretical calculations based on the trace 
formula derived by Balian and Duplantier. As there the 
polarization of the electric field is taken into account for 
each unstable periodic orbit individually, this good agree- 
ment is not in contradiction to the discrepancy obtained 
for the spectral properties. 




FIG. 5: (Color online) Upper panel: Nearest neighbor spacing 
distribution (histogram) of the billiard for ri = \/2?"2 com- 
pared to the GOE prediction (dashed line) . Lower panel: E 2 - 
statistics (thin line) compared to the GOE prediction (dashed 
line). 



metric under the symmetry operation and belong to dif- 
ferent irreducible representations (IR). Assuming that 
the billiard is chaotic, the spectral properties of levels 
within each IR are expected to coincide with those of 
random matrices from the GOE. However, the spectral 
statistics for the whole set of eigenvalues should coincide 
with RMT for a superposition of two independent GOEs. 
Let Hi and Hi be two random matrices drawn from 
the GOE and consider an ensemble of random matrices 
of the form [H [33|, [il 



H 



H x 
H 2 



(10) 



A. Symmetry breaking 

For r\ — r2, the billiard is symmetric with respect to 
reflections at the z = plane and a subsequent 7r/2- 
rotation about the z-axis. Accordingly, in this case the 
wave functions of the billiard are symmetric or antisym- 



For a generic chaotic system with two symmetry classes, 
the spectral properties of the eigenvalues associated with 
each symmetry class are given by those of the random 
matrices H\ and i? 2 , respectively, whereas the complete 
set of eigenvalues is described by those of H itself. Since, 
in our case, the number of eigenvalues associated with 



the two symmetry classes are (approximately) equal, the 
matrices Hi and H 2 are chosen of equal dimension for 
the theoretical description of the spectral properties of 
the quantum billiard. In Fig. [5] we compare the E 2 - 
statistics for the eigenvalues of the quantum billiard (thin 
line) with that of an ensemble of random matrices of the 
form Eq. (jTUJ) (dashed line), which is known analytically 
[34| . We observe significant deviations, which cannot be 
explained by an additional family of nongeneric modes. 
Indeed, Fig. [7] shows the fluctuating part of the stair- 
case function iVfl uc (/) for the case r\ — r 2 (thin line) 
and the contributions of the nongeneric modes resulting 
from Eq. ||7J) (thick line). As for the case r\/r2 = y2, 
the smooth oscillations of Na uc (f) are well described by 
our expression. From this we may conclude that the adi- 
abatic method described in Section |TT] yields a good ap- 
proximation for the contribution of the nongeneric modes 
to the staircase function also for r\ = r 2 . We also verified, 
that the deviations are not due to insufficient numerical 
accuracy in the computation of the eigenvalues. How can 
this puzzle be understood? 




FIG. 6: (Color online) £ 2 -statistics for the billiard with 
n — T2 (thin line) compared to the GOE prediction (dashed 
line) and to a random matrix model which also includes a 
Poissonian sequence of random levels (dots) . 

Recall that the billiard has a stable orbit for r\ = r 2 
and several marginally stable orbits besides the bouncing 
ball orbits. This suggests that their presence causes the 
deviation depicted in Fig. [6] In order to test this hypoth- 
esis, we considered the following random matrix ensem- 
ble. Each random matrix consists of two block matrices 
Hi and H 2 of equal dimension that are drawn from the 
GOE. These model the two symmetry classes of eigen- 
states associated with the chaotic part of the billiard 
(see Eq. pH|) h An additional diagonal random matrix 
of much smaller dimension models eigenstates associated 
with the few stable orbits. This diagonal matrix thus 
exhibits Poisson statistics. The dimension of the latter 
matrix was chosen equal to 25, that of Hi and H 2 equal 
to 250. The E 2 -statistics of this random matrix model 




FIG. 7: (Color online) The fluctuating part Nfi uc of the stair- 
case function for the billiard with n = r-2 (thin line) compared 
to the contributions from the nongeneric modes (thick line). 



is plotted as a dotted line in Fig [6] It agrees very well 
with that of the billiard. This suggests that the stable 
and marginally stable orbits are responsible for the ob- 
served deviations from the prediction of standard RMT. 
Most interestingly, the stable orbits disappear immedi- 
ately when the ratio ri/r 2 differs slightly from one. 

For geometries n ^ r%, the symmetry of the billiard is 
broken. The underlying quantum system "sees" this sym- 
metry breaking once the perturbation due to the symme- 
try breaking is of the order of the mean level spacing, or, 
alternatively, once the geometric distortions due to the 
symmetry breaking are of the scale of one wave length. 
For a completely broken symmetry, RMT predicts agree- 
ment of the spectral properties of the billiard with those 
of one GOE. Quantum mechanically, however, we found 
that for j*i ~ r 2 , the spectral properties coincide with 
those of random matrices applicable to chaotic systems 
with a partially broken symmetry [lj, [33|, [3J, [33, (36| , 



H = 



Hi 

H 2 



VXD 



V 
V T 



(11) 



In this random matrix model, the first (second) matrix 
preserves (breaks) the symmetry. The size of the sym- 
metry breaking is set by the dimensionless parameter A 
measured on the scale of the mean spacing D. Here as 
in Eq. (JTUJ) , Hi and H 2 are random matrices drawn from 
the GOE. The symmetry breaking is modeled by the off- 
diagonal blocks V and V T , where the random matrix V is 
real with no symmetries. For this random matrix model, 
the E 2 -statistics is known analytically for arbitrary val- 
ues of A. 

Figure [H] shows the E 2 -statistics (thin line) for an in- 
creasing ratio ri/r 2 of the radii of the billiard. For 
comparison, we also show the E 2 -statistics for one GOE 
(dashed line), for two GOEs (dotted line), and for the 
random matrix model (|11[) (thick line). The values of A 



given in the figure are obtained by a fit of the model (fTTj) 
to the data. For r 1 /r 2 = 1.0025 there is perfect agree- 
ment with the random matrix model Eq. (fTQ|) that de- 
scribes chaotic systems with a conserved symmetry. On 
the one hand, the ratio r\jr% — 1.0025 deviates (suffi- 
ciently strong) from one and the stable island has disap- 
peared. On the other hand, this ratio is still so close to 
one that the quantum mechanics is unable to resolve the 
symmetry breaking. For increasing values of the ratio 
r\jr 2l the symmetry breaking is revealed in the spectral 
statistics, and the parameter A in Eq. (Tl"Tj) increases from 
zero. Eventually, for r\jr 2 ~ 1.1 . . . 1.2, the E 2 -statistics 
approaches that of one GOE. This result is in agreement 
with a semiclassical estimate. The symmetry breaking 

We have 



3 

1 



< 



n. - r2 



is resolved for wave lengths 2ir/k 
maximal wave momenta kr 2 ~ 35 and can thus resolve 
symmetry breaking of the order r\fr 2 > 1 + 2ir/(kr 2 ) « 
1.18. Note that we can also base our semiclassical es- 
timate on the smooth part N amoot \ 1 (k) of the staircase 
function. This function is quadratic in the symmetry- 
breaking parameter (ri — r 2 ). At maximal wave momenta 
kr 2 rs 35 nonzero values of this parameter lead to signif- 
icant changes of Nn uc (k) for r\jr 2 > 1.18. In conclusion, 
the classically abrupt change of the symmetry properties 
of the system is accompanied by a gradual change of the 
spectral properties of the corresponding quantum system 
from those of chaotic systems consisting of a superposi- 
tion of two symmetry classes to those of chaotic systems 
with no further symmetries. 



We add here some more notes concerning the experi- 
ment with a microwave resonator [22j. There, the spec- 
tral properties were also described with the model given 
in Eq. (fTT|). where one of the block matrices Hi, H 2 de- 
picts the properties of the TE, the other those of the 
TM modes and deviations from GOE-behavior were in- 
terpreted as due to a partial decoupling of them. In the 
experiment the ratio r\/r 2 ~ \f2 with n = 200.0 mm 
and r 2 = 141.4 mm was kept fixed, while the value of 
the resonance frequency /, that is of k — =^J- with c the 
velocity of light, was varied up to / = 20 GHz. 



Is for this choice of the radii and of the frequency range 
the symmetry breaking discussed above observable? The 
breaking of the symmetry existent for T\ — r 2 is resolv- 
able for wavelengths smaller than \n — r 2 \, that is for 
27r/fc = c/ f < |t*i — 7*2! = 58.6 mm. Accordingly, for 
excitation frequencies / > 5.12 GHz good agreement of 
the spectral properties with those of random matrices 
from the GOE are expected. In the experiment, how- 
ever, deviations from one GOE were observed up to ap- 
proximately 17 GHz. Hence, they cannot be explained 
with the particular mechanism of symmetry breaking dis- 
cussed above. In order to resolve the discrepancy between 
theory and experiment numerical computations of the full 
vectorial Helmholtz equation are desired. 



— ' — 1 — ' — 1 — ' — r 

■r/r 2 =1.0000, XsO.OO 




FIG. 8: (Color online) E 2 -statistics for various ratios ri/ra of 
the billiard (thin line) compared to that of two GOEs (dotted 
line), one GOE (dashed line), and for the random matrix 
model in Eq. ([TO]) (thick line). 



B. High- lying states 

Strictly speaking, the semiclassical approximation and 
the Bohigas-Giannoni-Schmit [12, [37], [38|] conjecture, 
which states, that the spectral fluctuation properties in 
quantum systems with a chaotic classical dynamics co- 
incide with those of random matrices drawn from the 
Gaussian random-matrix ensembles, apply in the regime 
when the wave length is the smallest length scale of the 
system, i.e. kr 2 ^> 1 must hold. Nevertheless, we saw in 
the previous section that the semiclassical analysis of the 
length spectrum is also useful for the low- lying states. In 
this subsection, we compute high-lying quantum states 
of the billiard and perform level statistics. For the com- 
putation of high-lying levels, we employ Prosen's gener- 
alization [17J of the two-dimensional method by Vergini 
and Saraceno 3$] to three-dimensional billiards. This 
method is particularly suited for high-lying eigenstates 
and determines a stretch of eigenstates around some ar- 
bitrary wave momentum fco with kor\ 3> 1. This method 
yields 1669 levels in the regime 79 < kr 2 < 86. 

Figure [9] compares the spectral properties of the eigen- 
values in the lowest part of the spectrum (thin line) with 
those of the high-lying eigenvalues (dashed line) and the 
GOE (dashed-dotted line). While the E 2 -statistics of 
the low-lying set of eigenvalues agrees well with that of 
random matrices from the GOE, we observe significant 
deviations for that of the high-lying ones. This points to 
inadequacies of our procedure to extract the contribution 
of the nongeneric modes to the fluctuating part of the res- 
onance density. Recall that our semiclassical formula ([7]) 
takes account of the leading nongeneric modes. These are 
contributions that scale as k 3 with increasing wave mo- 
mentum. Next-to- leading order contributions scale as k 2 . 
These are only partly included. In our ansatz ([5]) we inte- 
grate over rectangles with a z-dependent area. There are 
many billiards of different shape (but identical volume) 



that have rectangular cross sections. However, evidently, 
our procedure is well suited for small or moderate values 
of kr-i , but it seems to be insufficient in the semiclassical 
regime. 




FIG. 9: (Color online) E 2 -statistics of the high-lying stretch 
of levels around kr^ ~ 80 (dashed-dotted line) compared to 
that of the low-lying levels (thin line) and the GOE prediction 
(dashed line). 



IV. SUMMARY 



We studied the classical and quantum mechanics of 
the three-dimensional stadium billiard. This billiard con- 
sists of two quarter cylinders that are rotated by 90 de- 
grees with respect to each other and is classically chaotic 
for unequal radii of the quarter cylinders. We studied 
the nongeneric modes of the billiard that are due to 
bouncing-ball orbits and orbits within a rectangular cross 
section of the billiard, and gave analytical expressions for 
their contribution to the staircase function. The quan- 
tum mechanical length spectrum can be understood in 
terms of the nongeneric modes and unstable periodic or- 
bits. For unequal radii of the two quarter cylinders, the 
level statistics agrees with predictions from random ma- 
trix theory. For equal radii, we found deviations from 
standard random matrix theory and these can be at- 
tributed to a stable and a few marginally stable orbits. 
We modeled a random matrix ensemble on this basis and 
found very good agreement with our data. We studied 
the level statistics as a function of the ratio of the two 
radii. For slightly unequal radii, there is a sudden tran- 
sition from a system with mixed phase space to a sys- 
tem with a chaotic phase space. For larger differences 
of the radii, we follow the symmetry-breaking transition 
towards the spectral statistics of one GOE. 



APPENDIX. FINITE ELEMENT METHOD 

In order to calculate a complete sequence of eigenvalues 
for the three-dimensional stadium billiard, we use a finite 
element method (see for example [40J], ch. 11). This 
discretization of the eigenvalue problem 



ipi 



kfipi in Cl, 
on dn, 



leads to a generalized eigenvalue problem 



A h $ = ^M h if>>, 



(12) 



(13) 



where A h ,M h G 



are finite dimensional matrices. 



The eigenvalues X\ and eigenfunctions tp^ of the dis- 
crete problem approximate the eigenvalues/functions of 
the continuous problem (fT2"|) . Here the parameter h G K. + 
represents the fineness of the test functions used in the 
discretization process. With the usual assumptions (see 
pOj). we get the error-estimate 



A? - Xi < CXl +1 h 



2/ 



(14) 



where I is the polynomial degree of the used test func- 
tions. This estimate shows that it is possible to compute 
an arbitrary number of eigenvalues with a given accuracy 
by using a sufficiently fine discretization. Using standard 
piecewise linear functions, we obtain convergence of order 
h 2 so that small values of h are required for high accuracy. 
Accordingly, the dimension of matrices A h and M h be- 
comes very large, and solving the generalized eigenvalue 
problem (|13|) is extremely costly. 

Using test functions of higher degree yields a signifi- 
cantly improved performance of the method. We suggest 
to use a special modification of tensor product B-splines 
with arbitrary coordinate degree. While, independent 
of I, B-splines may do with only one coefficient per el- 
ement, standard approaches of higher order require a 
much larger number, which is increasing with the degree. 
Since computation time to solve (|13[) is growing rapidly 
with the number of degrees of freedom, this aspect is 
crucial for achieving high accuracy. 

So far, despite these advantages, B-splines have rarely 
been used in finite element applications, the main rea- 
sons being Dirichlet boundary conditions and stability. 
A solution to the se p roblems was proposed by Hollig, 
Reif and Wipper [4JJ. Stability is achieved by linking 
potentially unstable B-splines near the boundary of the 
domain in an appropriate way to B-splines in the in- 
terior of the domain. Following old, but little known 
ideas of Kantorovitch and Krylow, homogeneous Dirich- 
let boundary conditions are ensured by multiplying all 
test functions by a weight function w which vanishes on 
the boundary and is positive inside the domain Q. It can 
be shown that the resulting weighted extended B-splincs 
(web-splines) form a stable basis and possess the same 
approximation order as the underlying tensor product 
splines space (il lil |43J . 
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In this work we apply the WEB-spline method to ap- 
proximate eigenvalues of the Laplacian on the three- 
dimensional stadium billiard O C I 3 . The implemen- 
tation includes a specifically designed high accuracy in- 
tegration algorithm. It is based on precomputed pro- 
jections of the boundary grid cells and iterated one- 
dimensional Gauss quadrature. 

We use tensor product B-splines of degree I = 4 on a 
grid with 25 x 25 x 40 cells. This leads to a generalized 
eigenvalue-problem (fl~3|) of dimension 9250. To deter- 
mine the eigenvalues, we use the Cholesky-decomposition 
M h = LL l to compute the matrix 



B h := L- 1 A h (L t ) 



(15) 



The eigenvalues of B are just the generalized eigenvalues 
of (A h ,M h ) so that standard software can be used to 
compute the Af. Compared with iterative methods, the 
advantage of this procedure is that one can be sure that 
no eigenvalues in the domain of interest are lost. On the 
other hand, the maximal dimension that can be handled 
that way is limited by the fact that B h is a dense matrix. 
Of course, only a certain fraction of the eigenvalues A^ 
of B h provides a reasonable approximation of some Aj. 
The estimate (fT4"|) suggests that smaller eigenvalues are 



approximated better than larger ones, but it does not 
provide actual bounds since the constant C is not known 
explicitly. Computations on other domains, where the 
exact eigenvalues A^ are explicitly known, indicate that 
at least the first 16% of the eigenvalues Aj 1 , when ordered 
my modulus, have a relative error less than 0.1% For 
the statistical evaluation in this paper, 1200 out of 9250 
eigenvalues, i.e., 13% are used. 
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